In [ ]:
import numpy as np  
import matplotlib.pyplot as plt  
import scipy 
import math  


# finite difference method to calculate the second derivative of a function
def d2f(f, dz, Nz):
    sol = np.zeros((Nz))  # Initialize solution array
    for i in range(1, Nz - 2):
        sol[i] = (f[i - 1] - 2 * f[i] + f[i + 1]) / (dz * dz)
    # Apply boundary conditions at i=0
    sol[0] = (2 * f[0] - 5 * f[1] + 4 * f[2] - f[3]) / (dz * dz)
    # Apply boundary conditions at i=Nz-1
    sol[Nz - 1] = (2 * f[Nz - 1] - 5 * f[Nz - 2] + 4 * f[Nz - 3] - f[Nz - 4]) / (dz * dz)
    return sol

def df(f,dz,Nz):
    sol = np.zeros((Nz))
    i=0
    sol[i]=((2*f[i])-(5*f[i+1])+(4*f[i+2])-f[i+3])/(dz*dz)
    i=Nz-1
    sol[i]=((2*f[i])-(5*f[i-1])+(4*f[i-2])-f[i-3])/(dz*dz)
    for i in range(1,Nz-2):
        sol[i]=(f[i-1]+(-2*f[i])+f[i+1])/(dz*dz)
    return sol

def solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D):
    #intials()
    Br=[]
    Bphi=[]
    decay=[]

    if 0 in time_plot:
        Br.append(Br_t0)
        Bphi.append(Bphi_t0)

    
    for j in range(0, Nt + 1):
        
        k1r = dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,dt*j)
        k1phi = dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,dt*j)
        
        k2r = dBrdt(Br_t0+k1r*dt/2, Bphi_t0+k1phi*dt/2, z,D,alpha,dz,Nz,vzfn,dt*j)
        k2phi = (dBphidt(Br_t0+k1r*dt/2, Bphi_t0+k1phi*dt/2, z,D,alpha,dz,Nz,vzfn,dt*j)) 
        
        k3r = (dBrdt(Br_t0+k2r*dt/2,Bphi_t0+k2phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        k3phi = (dBphidt(Br_t0+k2r*dt/2,Bphi_t0+k2phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        
        k4r = (dBrdt(Br_t0+k3r*dt,Bphi_t0+k3phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        k4phi = (dBphidt(Br_t0+k3r*dt/2,Bphi_t0+k3phi*dt,z,D,alpha,dz,Nz,vzfn,dt*j))
        Br_t0 = Br_t0 + ((dt / 6.0) * (k1r +( 2 * k2r) + (2 * k3r) + k4r))
        
        Bphi_t0 = Bphi_t0 + ((dt / 6.0) * (k1phi + (2 * k2phi) + (2 * k3phi) + k4phi))
               
        #(Anti Symmetric Ghost Zone)
        Bphi_t0[1]=Bphi_z0
        Br_t0[Nz-2]=Br_zf
        Br_t0[1]=Br_z0
        Bphi_t0[Nz-2]=Bphi_zf
        Br_t0[Nz-1]=Br_z0-Br_t0[Nz-3]
        Bphi_t0[Nz-1]=Bphi_z0-Bphi_t0[Nz-3]
        Br_t0[0]=Br_z0-Br_t0[2]
        Bphi_t0[0]=Bphi_z0-Bphi_t0[2]

        decay.append(np.log10(np.sqrt((np.mean(Br_t0))**2+(np.mean(Bphi_t0))**2)))

        if j in time_plot:
            Br.append(Br_t0)
            Bphi.append(Bphi_t0)
            #B_pitch.append(np.copy(B))

    pitch=np.arctan(np.array(Br)/np.array(Bphi))
    return Br,Bphi,pitch,decay


def plotfn(Br, Bphi, pitch, decay, time_plot,t,z,Nz):
    m, b = np.polyfit(t[2000:], decay[2000:], 1)

    fig, ax = plt.subplots(2, 2, figsize=(15, 10))  
    
    for i in range(len(time_plot)):
        ax[0, 0].plot((z[1:Nz-1])/5, Br[i][1:Nz-1], label=f't={time_plot[i]}')
        ax[0, 1].plot((z[1:Nz-1])/5, Bphi[i][1:Nz-1], label=f't={time_plot[i]}')
        ax[1, 1].plot((z[1:Nz-1])/5, pitch[i][1:Nz-1], label=f't={time_plot[i]}')
    
    ax[0, 0].set_xlabel(r'z/$h_0$')
    ax[0, 0].set_ylabel(r'$B_r/B_0$')
    ax[0, 0].set_title(r'$B_r/B_0$ vs z')
    ax[0, 0].legend()
    
    ax[0, 1].set_xlabel(r'z/$h_0$')
    ax[0, 1].set_ylabel(r'$B_{\phi}/B_0$')
    ax[0, 1].set_title(r'$B_{\phi}/B_0$ vs z')
    ax[0, 1].legend()

    ax[1, 0].scatter(t, decay)
    ax[1, 0].plot(t[2000:], (m*t+b)[2000:],color='Red', label='linear fit')
    ax[1, 0].set_xlabel('$t/t_0$')
    ax[1, 0].set_ylabel(r'$B/B_0$')
    #ax[1, 0].set_title('exponential decay rate')
    
    #ax[1, 1].plot(z, pitch)
    ax[1, 1].set_xlabel(r'z/$h_0$')
    ax[1, 1].set_ylabel('pitch angle')
    ax[1, 1].set_title('pitch angle vs z')
    
    plt.tight_layout()
    filename = f'plot_{vzfn(5)}.png'  # Constructing the filename
    plt.savefig(f'images/{filename}')  # Saving the plot with the constructed filename
    print(f'Plot saved as {filename}')
    plt.show()
    plt.clf()
    # print("The Decay rate is",np.round((m*np.log(10)),4))
    print("The Decay rate is",np.round((m),4))

constant Vz = 0.80¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))


# Defning the velocity function
def vzfn(tv):
    return 0.80


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))




def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 50 
Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_32072\2457440500.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))
No description has been provided for this image 
The Decay rate is 0.0055

<Figure size 640x480 with 0 Axes>

Vz=0¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))


# def vzfn(tv):
#     return tv/(4.1*np.sqrt(tv+0.001))

# def vzfn(tv):
#     return np.sin(tv/80)
def vzfn(tv):
    return 0


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 50 
Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_27240\2703013185.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))
No description has been provided for this image 
The Decay rate is -0.0367

<Figure size 640x480 with 0 Axes>

$V_z$ is a linear function of time¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))



def vzfn(tv):
    return tv/51


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 50 
Nz = 21  # Number of spatial grid points
Nt = 15000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_32072\2457440500.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))
No description has been provided for this image 
The Decay rate is -0.0057

<Figure size 640x480 with 0 Axes>

$V_z(t) = \frac{L}{1 + e^{-k(t - t_0)}}$¶

In [ ]:
# def logistic_function(t, L, k, t0):
#     return L / (1 + np.exp(-k * (t - t0)))

# # Define the range of t values
# t_values = np.linspace(0, 50, 100)

# # Choose parameters for the logistic function
# L = 50  # Maximum value
# k = 0.1  # Controls the steepness of the curve
# t0 = 0  # Time of maximum growth

# # Compute the corresponding function values
# function_values = logistic_function(t_values, L, k, t0)


def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))


# def vzfn(tv):
#     return tv/(4.1*np.sqrt(tv+0.001))

# def vzfn(tv):
#     return np.sin(tv/80)
L = 0.874 # Maximum value
k = 0.1  # Controls the steepness of the curve
t0 = 0  # Time of maximum growth
def vzfn(tv):
    return L / (1 + np.exp(-k * (tv - t0)))


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 70 
Nz = 21  # Number of spatial grid points
Nt = 15000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_32072\2457440500.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))
No description has been provided for this image 
The Decay rate is 0.0043

<Figure size 640x480 with 0 Axes>

for D=-5 with same function as above¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))


# def vzfn(tv):
#     return tv/(4.1*np.sqrt(tv+0.001))

# def vzfn(tv):
#     return np.sin(tv/80)
L = 0.874 # Maximum value
k = 0.1  # Controls the steepness of the curve
t0 = 0  # Time of maximum growth
def vzfn(tv):
    return L / (1 + np.exp(-k * (tv - t0)))


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 70 
Nz = 21  # Number of spatial grid points
Nt = 15000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# Defining D and alpha
D, alpha = -5, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_27240\1860986583.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_0.5440294554704209.png
No description has been provided for this image 
The Decay rate is 0.0042

<Figure size 640x480 with 0 Axes>

for D=-50¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))



L = 0.874 # Maximum value
k = 0.1  # Controls the steepness of the curve
t0 = 0  # Time of maximum growth
def vzfn(tv):
    return L / (1 + np.exp(-k * (tv - t0)))


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))




def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))



z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 70 
Nz = 21  # Number of spatial grid points
Nt = 15000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



#definig D and alpha
D, alpha = -50, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_27240\1860986583.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_0.5440294554704209.png
No description has been provided for this image 
The Decay rate is 0.0043

<Figure size 640x480 with 0 Axes>

When $V_z$ is large

$V_z$=t/37

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))




def vzfn(tv):
    return tv/37

def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 42 
Nz = 21  # Number of spatial grid points
Nt = 100000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_27240\1860986583.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_0.13513513513513514.png
No description has been provided for this image 
The Decay rate is 0.0054

<Figure size 640x480 with 0 Axes>

when $V_z$ is not large enough¶

In [ ]:
def vzfn(tv):
    return 0.5

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_32072\2457440500.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))
No description has been provided for this image 
The Decay rate is -0.0031

<Figure size 640x480 with 0 Axes>

$V_z$ is a function of z¶

$V_z$= 0.8/5*z

In [ ]:
def plotfn(Br, Bphi, pitch, decay, time_plot,t,z,Nz):
    m, b = np.polyfit(t[2000:], decay[2000:], 1)

    fig, ax = plt.subplots(2, 2, figsize=(15, 10))  
    
    for i in range(len(time_plot)):
        ax[0, 0].plot((z[1:Nz-1])/5, Br[i][1:Nz-1], label=f't={time_plot[i]}')
        ax[0, 1].plot((z[1:Nz-1])/5, Bphi[i][1:Nz-1], label=f't={time_plot[i]}')
        ax[1, 1].plot((z[1:Nz-1])/5, pitch[i][1:Nz-1], label=f't={time_plot[i]}')
    
    ax[0, 0].set_xlabel(r'z/$h_0$')
    ax[0, 0].set_ylabel(r'$B_r/B_0$')
    ax[0, 0].set_title(r'$B_r/B_0$ vs z')
    ax[0, 0].legend()
    
    ax[0, 1].set_xlabel(r'z/$h_0$')
    ax[0, 1].set_ylabel(r'$B_{\phi}/B_0$')
    ax[0, 1].set_title(r'$B_{\phi}/B_0$ vs z')
    ax[0, 1].legend()

    ax[1, 0].scatter(t, decay)
    ax[1, 0].plot(t[2000:], (m*t+b)[2000:],color='Red', label='linear fit')
    ax[1, 0].set_xlabel('$t/t_0$')
    ax[1, 0].set_ylabel(r'$B/B_0$')
    #ax[1, 0].set_title('exponential decay rate')
    
    #ax[1, 1].plot(z, pitch)
    ax[1, 1].set_xlabel(r'z/$h_0$')
    ax[1, 1].set_ylabel('pitch angle')
    ax[1, 1].set_title('pitch angle vs z')
    
    plt.tight_layout()
    filename = f'plot_{vzfn(5,2)}.png'  # Constructing the filename
    print(filename)
    plt.savefig(f'images/{filename}')  # Saving the plot with the constructed filename

    plt.show()
    plt.clf()
    # print("The Decay rate is",np.round((m*np.log(10)),4))
    print("The Decay rate is",np.round((m),4))




def Bphi_t0(z):
    return np.cos((np.pi/2)*z/5)

def Br_t0(z):
    return np.cos((np.pi/2)*z/5+(np.pi))




vz=0.8

def vzfn(tv,z):
    return vz*z/5


def alphafn(z):
    return (2*np.sin(np.pi*z/1000000))

# def alphafn(z):
#     return (np.sin(np.pi*z/0.761))


def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv,z)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv,z)*Bphi_t0,dz,Nz))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
z0 = -5  # start of spatial region in z
zf = 5  # end of spatial domain in z

t0=0.0
tf = 50 
Nz = 21  # Number of spatial grid points
Nt = 20000  # Number of time steps
dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size



z = np.linspace(z0, zf, Nz)
t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]



# D, alpha = -5, 1
D, alpha = -17.10, 1

plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_27240\3479477426.py:73: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

plot_0.32.png
No description has been provided for this image 
The Decay rate is -0.0285

<Figure size 640x480 with 0 Axes>
In [ ]:

updated with sixth orfer finite difference¶

In [ ]:
import numpy as np  
import matplotlib.pyplot as plt  
import scipy 
import math  
from tqdm import tqdm


# 6th order finite difference method to calculate the second derivative of a function
def d2f(f, dz, Nz):
    
    sol = np.zeros(Nz)  # Initialize solution array
    for i in range(3, Nz - 4):

        sol[i] = ((2 * f[i - 3]) - (27 * f[i - 2]) + (270 * f[i - 1]) - (490 * f[i]) + (270 * f[i + 1]) - (27 * f[i + 2]) + (2 * f[i + 3])) / (180 * (dz ** 2))

    

    return sol

def df(f,dz,Nz):
    # np.copy(f)
    sol = np.zeros(Nz)
    for i in range(3,Nz-4):
        #sol[i] = ((-2 * f[i - 5]) + (25 * f[i - 4]) - (150 * f[i - 3]) + (600 * f[i - 2]) - (2100 * f[i - 1]) + (2100 * f[i + 1]) - (600 * f[i + 2]) + (150 * f[i + 3]) - (25 * f[i + 4]) + (2 * f[i + 5])) / (2520 * dz)
        sol[i] = ((- f[i - 3]) + (9 * f[i - 2]) - (45 * f[i - 1]) + (45 * f[i + 1]) - (9 * f[i + 2]) + (f[i + 3])) / (60 * dz)
    #sol=sol[3:-3]   
    return sol

def solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D):
    #intials()
    Br=[]
    Bphi=[]
    decay=[]

    if 0 in time_plot:
        Br.append(Br_t0)
        Bphi.append(Bphi_t0)

    
    for j in range(0, Nt + 1):
        
        k1r = dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,dt*j)
        k1phi = dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,dt*j)
        
        k2r = dBrdt(Br_t0+k1r*dt/2, Bphi_t0+k1phi*dt/2, z,D,alpha,dz,Nz,vzfn,dt*j)
        k2phi = (dBphidt(Br_t0+k1r*dt/2, Bphi_t0+k1phi*dt/2, z,D,alpha,dz,Nz,vzfn,dt*j)) 
        
        k3r = (dBrdt(Br_t0+k2r*dt/2,Bphi_t0+k2phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        k3phi = (dBphidt(Br_t0+k2r*dt/2,Bphi_t0+k2phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        
        k4r = (dBrdt(Br_t0+k3r*dt,Bphi_t0+k3phi*dt/2,z,D,alpha,dz,Nz,vzfn,dt*j)) 
        k4phi = (dBphidt(Br_t0+k3r*dt/2,Bphi_t0+k3phi*dt,z,D,alpha,dz,Nz,vzfn,dt*j))
        Br_t0 = Br_t0 + ((dt / 6.0) * (k1r +( 2 * k2r) + (2 * k3r) + k4r))
        
        Bphi_t0 = Bphi_t0 + ((dt / 6.0) * (k1phi + (2 * k2phi) + (2 * k3phi) + k4phi))
               
        #(Anti Symmetric Ghost Zone)
        Bphi_t0[3]=Bphi_z0
        Br_t0[3]=Br_z0

        Bphi_t0[Nz-4]=Bphi_zf
        Br_t0[Nz-4]=Br_zf

        Br_t0[Nz-1]=Br_zf-Br_t0[Nz-7]
        Bphi_t0[Nz-1]=Bphi_zf-Bphi_t0[Nz-7]

        Br_t0[Nz-2]=Br_zf-Br_t0[Nz-6]
        Bphi_t0[Nz-2]=Bphi_zf-Bphi_t0[Nz-6]

        Br_t0[Nz-3]=Br_zf-Br_t0[Nz-5]
        Bphi_t0[Nz-3]=Bphi_zf-Bphi_t0[Nz-5]

        Br_t0[0]=Br_z0-Br_t0[6]
        Bphi_t0[0]=Bphi_z0-Bphi_t0[6] 

        Br_t0[1]=Br_z0-Br_t0[5]
        Bphi_t0[1]=Bphi_z0-Bphi_t0[5]

        Br_t0[2]=Br_z0-Br_t0[4]
        Bphi_t0[2]=Bphi_z0-Bphi_t0[4]


        


          



        decay.append(np.log10(np.sqrt((np.mean(Br_t0))**2+(np.mean(Bphi_t0))**2)))

        if j in time_plot:
            Br.append(Br_t0)
            Bphi.append(Bphi_t0)
            #B_pitch.append(np.copy(B))

    pitch=np.arctan(np.array(Br)/np.array(Bphi))
    return Br,Bphi,pitch,decay


def plotfn(Br, Bphi, pitch, decay, time_plot,t,z,Nz):
    m, b = np.polyfit(t[2000:], decay[2000:], 1)

    fig, ax = plt.subplots(2, 2, figsize=(15, 10))  
    
    for i in range(len(time_plot)):

        ax[0, 1].plot(z[3:Nz-3], Bphi[i][3:Nz-3], label=f't={time_plot[i]/Nt*tf}')
        ax[1, 1].plot(z[3:Nz-3], pitch[i][3:Nz-3], label=f't={time_plot[i]/Nt*tf}')
        ax[0, 0].plot(z[3:Nz-3], Br[i][3:Nz-3], label=f't={time_plot[i]/Nt*tf}')
    
    ax[0, 0].set_xlabel(r'$z/h_0$')
    ax[0, 0].set_ylabel(r'$Br/B_0$')
    ax[0, 0].set_title(r'$Br/B_0 vs z/h_0$')
    ax[0, 0].legend()
    
    ax[0, 1].set_xlabel(r'$z/h_0$')
    ax[0, 1].set_ylabel(r'$B_{\phi}/B_0$')
    ax[0, 1].set_title(r'$B_{\phi}/B_0$ vs $z/h_0$')
    ax[0, 1].legend()

    ax[1, 0].scatter(t, decay)
    ax[1, 0].plot(t[2000:], (m*t+b)[2000:],color='Red', label='linear fit')
    ax[1, 0].set_xlabel(r'$t/t_0$')
    ax[1, 0].set_ylabel('log(B)')
    ax[1, 0].set_title('exponential decay rate')
    
    #ax[1, 1].plot(z/h_0$, pitch)
    ax[1, 1].set_xlabel(r'$z/h_0$')
    ax[1, 1].set_ylabel('pitch angle')
    ax[1, 1].set_title(r'$pitch angle vs z/h_0$')
    
    plt.tight_layout()
    filename = f'plot_{vzfn(5)}.png'  # Constructing the filename
    plt.savefig(f'images/{filename}')  # Saving the plot with the constructed filename
    print(f'Plot saved as {filename}')
    plt.show()
    print("The Decay rate is",np.round((m*np.log(10)),4))



def D_c(solve2,precision):
    D=0
    h=2
    while h>precision:  
        drate=-1 #dummy value
        print('current |D| =',D)
        while drate<0:
            Br,Bphi,pitch,decay=solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,-D)
            m, c = np.polyfit(t, decay , 1)
            drate=m*np.log(10)
            print('|D| =',D,' | decay rate =',drate)
            D=D+h   
        D=D-2*h
        h=h/2
    print("Critical dynamo number, |D_c| =",D+h/2)
    #return D+h/2

without $V_z$ outflow¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/h)

def Br_t0(z):
    return np.cos((np.pi/2)*z/h+(np.pi))


def vzfn(tv):
    return 0.01




def alphafn(z):
    return (np.sin(np.pi*z/h))

def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0))


# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
# z0 = -0.768 # start of spatial region in z
# zf = 0.768 # end of spatial domain in z

z0 = -1.1 # start of spatial region in z
zf = 1.1 # end of spatial domain in z

t0=0.0
tf = 2 

Nz = 51  # Number of spatial grid points
Nt = 10000  # Number of time steps

dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size
h=0.768


z = np.linspace(z0, zf, Nz)

t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]
#time_plot=[10000]



D, alpha = -20, 1
plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_1948\2908489377.py:97: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_0.01.png
No description has been provided for this image 
The Decay rate is 0.8629

$V_z$ = 1.045¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/h)

def Br_t0(z):
    return np.cos((np.pi/2)*z/h+(np.pi))


# def vzfn(tv):
#     return 0




def alphafn(z):
    return (np.sin(np.pi*z/h))


def vzfn(tv):
    return 1.045



def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))

# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
# z0 = -0.768 # start of spatial region in z
# zf = 0.768 # end of spatial domain in z

z0 = -1.1 # start of spatial region in z
zf = 1.1 # end of spatial domain in z

t0=0.0
tf = 2 

Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps

dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size
h=0.768


z = np.linspace(z0, zf, Nz)

t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[1,int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]
#time_plot=[10000]



D, alpha = -20, 1
plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_1948\2908489377.py:97: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_1.045.png
No description has been provided for this image 
The Decay rate is 0.0022

$V_z$ is a function of time¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/h)

def Br_t0(z):
    return np.cos((np.pi/2)*z/h+(np.pi))


# def vzfn(tv):
#     return 0




def alphafn(z):
    return (np.sin(np.pi*z/h))

L = 1.045 # Maximum value
k = 4  # Controls the steepness of the curve
t0 = 0  # Time of maximum growth

def vzfn(tv):
    return L / (1 + np.exp(-k * (tv - t0)))

# def vzfn(tv):
#     return 1.045



def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))



z0 = -1.1 # start of spatial region in z
zf = 1.1 # end of spatial domain in z

t0=0.0
tf = 2 

Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps

dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size
h=0.768


z = np.linspace(z0, zf, Nz)

t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[1,int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]
#time_plot=[10000]



D, alpha = -20, 1
plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_1948\2908489377.py:97: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_1.0449999978460944.png
No description has been provided for this image 
The Decay rate is 0.011

For $V_z$ = 2¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/h)

def Br_t0(z):
    return np.cos((np.pi/2)*z/h+(np.pi))


# def vzfn(tv):
#     return 0




def alphafn(z):
    return (np.sin(np.pi*z/h))


def vzfn(tv):
    return 2



def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))

# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
# z0 = -0.768 # start of spatial region in z
# zf = 0.768 # end of spatial domain in z

z0 = -1.1 # start of spatial region in z
zf = 1.1 # end of spatial domain in z

t0=0.0
tf = 2 

Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps

dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size
h=0.768


z = np.linspace(z0, zf, Nz)

t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]
#time_plot=[10000]



D, alpha = -20, 1
plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_1948\2908489377.py:97: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_2.png
No description has been provided for this image 
The Decay rate is -0.6786

$V_z$ =0.5¶

In [ ]:
def Bphi_t0(z):
    return np.cos((np.pi/2)*z/h)

def Br_t0(z):
    return np.cos((np.pi/2)*z/h+(np.pi))


# def vzfn(tv):
#     return 0




def alphafn(z):
    return (np.sin(np.pi*z/h))


def vzfn(tv):
    return 0.5



def dBrdt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Br_t0,dz,Nz)-df(alpha*alphafn(z)*Bphi_t0,dz,Nz)-df(vzfn(tv)*Br_t0,dz,Nz))
def dBphidt(Br_t0,Bphi_t0,z,D,alpha,dz,Nz,vzfn,tv):
    return (d2f(Bphi_t0,dz,Nz)+((D/alpha)*Br_t0)-df(vzfn(tv)*Bphi_t0,dz,Nz))

# z0 = 0.0  # start of spatial region in z
# zf = 10.0  # end of spatial domain in z
# z0 = -0.768 # start of spatial region in z
# zf = 0.768 # end of spatial domain in z

z0 = -1.1 # start of spatial region in z
zf = 1.1 # end of spatial domain in z

t0=0.0
tf = 2 

Nz = 21  # Number of spatial grid points
Nt = 10000  # Number of time steps

dz = (zf - z0) / (Nz - 1)  # Spatial step size
dt = tf / Nt  # Time step size
h=0.768


z = np.linspace(z0, zf, Nz)

t = np.linspace(t0, tf, Nt+1)

Br_t0=Br_t0(z)
Bphi_t0=Bphi_t0(z)

Bphi_z0=0
Bphi_zf=0
Br_z0=0
Br_zf=0

#print(Bphi_z0,Bphi_zf,Br_z0,Br_zf)
time_plot=[int(Nt/10),int(Nt/5),int(Nt/2),int(Nt/1.5),int(Nt/1.2),int(Nt/1.1),Nt]
#time_plot=[10000]



D, alpha = -20, 1
plotfn(*solve2(Br_t0,Bphi_t0,Nz,Nt,time_plot,Bphi_z0,Bphi_zf,Br_z0,Br_zf,D),time_plot,t,z,Nz)

C:\Users\junaid\AppData\Local\Temp\ipykernel_1948\2908489377.py:97: RuntimeWarning: invalid value encountered in divide
  pitch=np.arctan(np.array(Br)/np.array(Bphi))

Plot saved as plot_0.5.png
No description has been provided for this image 
The Decay rate is 0.1993
In [ ]:
L = 1.045 # Maximum value
k = 4 # Controls the steepness of the curve
t0 = 0  # Time of maximum growth

def vzfn(tv):
    return L / (1 + np.exp(-k * (tv - t0)))
tv=np.linspace(0,2,100)
a=vzfn(tv)
plt.plot(tv,a)
Out[ ]:
[<matplotlib.lines.Line2D at 0x1f6852b6240>]
No description has been provided for this image